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Abstract: In This paper we illustrate new methods of online nonlinear estimation 
applied to the lateral deflection of an elastic beam from on board measurements 
of angular rates and angular accelerations. We contrast the development of the 
filter equations, together with practical issues of their numerical solution as devel- 
oped from global linearization by nonlinear output injection with the usual method 
of the the extended Kalman filter(EKF).We show how nonlinear estimation due 
to gyroscopic coupling can be implemented as an adptive covariance filter using 
off-the-shelf Kalman filter alghorithms. The effect of the global linearization by 
nonlinear output injection is to introduce a change of coordinates in which only 
the process noise covariance is to be updated in online implementation. This is in 
contrast to the computational approach which arises in EKF methods arising by 
local linearization with respect to the current conditional mean. We also highlight 
processing refinements for nonlinear estimation based on optimal ,nonlinear inter- 
polation between observations. In these methods the extrapolation of the process 
dynamics between measurement updates is obtained by replacing a transition ma- 
trix with an operator spline that is optimized off-line from responses to selected 
test inputs . 
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1 Introduction 


In this paper, a new nonlinear, nonparametric method is proposed for off- 
line modeling and on-line estimation of nonlinear dynamic systems. For 
illustration, it is applied to the estimation of the deformation of an elastic 
structure, undergoing rapid rotational maneuvers. 

In these circumstances, the structural stiffness and damping coefficients 
depend on the angular acceleration u,the angular rate u and the square of 
the angular rate .In the single axis case, the excitation of the structure is 
represented by the vector u T = (w, a? 2 , 2a;), to which the structural dynam- 
ics responds as a ’’bilinear” ( i.e., parametrically excited) system. A similar 
technique for multiaxial rotations yields a bilinear model with respect to 
matrix valued excitations. 

Two methods of estimation and modeling are combined to achieve de- 
formation state determination: 


• A method based on a feedback linearized procedure which gives an esti- 
mate by a filter constructed from the equivalent linear dynamics, which 
is faster than the extended Kalman filter. 

• The modeling of the deformation state of the structure by means of a 
Volterra series interpolator. 

2 Simplified Model of a Deformable Structure and Equations 
of Motion 

For purposes of illustration of the principles involved, the structure will con- 
sist of a primary mirror, attached to a spacecraft (containing the hardware 
of the slewing controller), and a secondary mirror attached to the central 
one in the shape of a Cassegrain telescope by means of massless links. The 
primary mirror structure will also be regarded as attached to the spacecraft 
by means of a massless link. Equivalently the same model can be thought to 
represent a laser beam expander, as in Figure 1 . More realistic models, such 
as in [1] or [2], exibit the same form of interaction between the rotational 
and vibrational dynamics. The simplified telescope part of the structure 
can itself be modeled as a system of two masses attached together by a sin- 
gle ’’equivalent” link, with ’’equivalent” stiffness and damping coefficients, so 
that the same restoring and dissipation forces at the secondary are obtained 
as if with more than one link. The modeling of such a deformable body is 



summarized in Figure 2 . 

If one takes now only the vibrational equation of motion ,and set the rota- 
tion around a single axis b 3 ,i.e, J =6b 3 ,and if the translational acceleration 
term is substituted from the translational equation, then one finds: 

My + (C + 26JM)y + (K + (0 + j6 2 )JM)y + (0 + 6 2 J)JM P = f (1) 

where J is not an inertia term, but rather an augmented symplectic matrix, 

made of blocks ^ ^ while M is a modified structure mass matrix to 

account for contributions due to translation, and p\s the 2nxl matrix compo- 
nents of the vectors pi from the undeformed appendage mass centers. Here 
y denotes the (2n x 1) (for planar motion or(3n X 1) for out of plane motion) 
matrix of deflection coordinates of the center of mass of n appendages from 
their undeformed positions, n = 2 in the case of the secondary mirror and 
the spacecraft platform regarded as appendages of the primary, 
f is the (2n x 1) (say, 4x1 here)matrix of body coordinates of external 
forces acting on the centers of mass of the n appendages. 

All other notations used are found in [12], [13]. 


Let now rj be a new variable such that 


T] = M(y + p) 

(2) 

and let 

C = CM- 1 

(3) 

K = KM- 1 

0) 

as well as 

u T = {u,u 2 ,2lo) 

(5) 

Then the vibrational equation of motion becomes: 


^ + ( C + U3«/)?7 + u x J + U3</ 2 )t7 = f + K p 

(6) 

This transformed equation can also be written in the bilinear form , which 
will be used frequently in the following sections, 

X =/lX + 5(X)u + b 

(7) 

where: 

A - ( to] m \ 

\{-K) [~C]J 

(8) 
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while Xjand X 2 are the vector components of the state vector X = {r^ ,r) T ) T . 

This simplified model,insofar as the links are regarded to be mass— less , 
exibits all the coupling effects between slewing motion and vibrational mo- 
tion. A distributed model under the assumption of symmetry about the mass 
center also yields product terms between 6 2 and structural deformations, and 
can be found in Chapter 9 of the book by Junkins and Turner[3]. That model 
of a symmetric four appendage spacecraft can also be used to illustrate the 
procedures being developed in this study, if desired, although damping must 
be present, so that the matrix A above will be stable. A slewing Timo- 
shenko beam model likewise exhibits this gyroscopic coupling effect, and 
also accounts for damping, so that A becomes stable, as found in [4]. 

3 Estimation of the state by means of observers 

This part of the estimation technique will deal mainly with updating the 
state from sensor data. 


3.1 Extended Kalman filter formulation 

Usually, when one deals with a nonlinear system of wdiich the state variables 
cannot all be observed (or are corrupted with noise) ,then the most com- 
monly used method of filtering or smoothing is the extended Kalman filter 
formulation [5]. Let the dynamical system be modeled as shown, 


x = f(x) + G(x)u* + G(x)£ 
y = h(x + v) 


where u* is the deterministic(mean) part of the input, £ is a zero mean input 
noise, and h is defined in our case as h = {h\ ,h^) T , where 


f M-) = A/ -1 (. - n) 

{ MO = 


( 12 ) 


Let Ri - r) = E[ui(t)vi(r) T ] be the covariance matrix of the sensor noise 
vectors 1 /,, with R = diag(Ri,R 2 ) and let Q6(t - r) = .E[£(i)f(r) T ] be 
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the covariance matrix of the actuator noise with V{ and £ presumed to be 
uncorrelated for simplicity. 

The propagation error matrix is defined by P , which satisfies the following 
Riccati differential equation along the estimated state trajectory x: 

P = [ jfk P + P\%\1 + 

'£Sr'l£gr 

x can be expressed as an observer, 

f x = f(x) + G(x)u* + K(y - y) 

\ y = K&) 

where K is the extended Kalman gain , and is defined as follows: 

K = p \^)l aix R- 1 (15) 

A procedure [6] based on a change of variables, in preliminary studies gave 
a faster computation time .This procedure is outlined next. 


(13) 

(14) 


3.2 Feedback linearized procedure: 


The idea is to change the state configuration of the original system , which 
has the particular form below: 


ii = Fi(x i)x 2 

x 2 = F 2 (x)(u* + £) + f 2 (x) 

(16) 

{ V\ =hi(xi + vi) 
l J /2 = h 2 (x 2 + V 2 ) 

(17) 


By using the change of variables x[ = x\,x 2 = F\(x\)x 2 , 

«' = *! = F 2 (x)u* + f 2 (x) and y[ = = K l {y 2 ) one gets 


( x 'i = x 'l 

\ x> 2 = u'* + ? 


(18) 


where = F\ F 2 £ and v[ = F \ ,so that the covariance of is approximated 

by 


Q' = F 1 (x 1 )F 2 (x)QF 2 (x) T F 1 (x 1 ) T 


(19) 



while that for i/[ is approximated by 


R'i = F l (i 1 )R i F 1 (x 1 ) T (20) 

Then the new error covariance matrix propagation is derived from the fol- 
lowing Riccati differential equation: 


P>-(V>) WW + p-fW [0]'l 
uo] m) F+F \\i\ io]J 

(j$j)<?'([o] W)- 
uo] w) uo] m) 


(21) 


The observed deformation state is also propagated in the usual manner: 

'»-<(B 18-UB !»* 

where the innovation process gain is now given as follows: 


A" = P'(Mo)W 1 (23) 

For the case of single axis slew-induced structural deformation estimation 
one has F\ = / and F 2 = lower half of B( X) defined by equation (2). In 
particular, one finds R\ = ,so that,in contrast to the extended Kalman 
filter, only the P'~ independent forcing term of the equation (21) given by 
Q * to be updated, all coefficients being now constants. In dealing with 
this procedure a 25% increase in speed, with comparable accuracy has been 
found in preliminary simulations discussed below. The problem in using 
either one of those two estimation techniques (even the faster one) in more 
realistic models than the example used here, is the high dimensionality of 
the filter then required, which may not be accommodated by the on-board 
data processing rate, causing estimation delays. 



3.3 Estimation Examples: 


A simplified model of a beam expander was represented by a primary mirror 
mass elastically linked to a secondary mirror mass according to the simpli- 
fied model outlined in Section 2. Restoring forces and dissipative forces 
proportional to relative secondary mirror motion were modeled at the sec- 
ondary.The structural parameters are found in the companion example in 
Subsection 5.1. A piecewise constant angular acceleration was commanded, 
representing the acceleration-deceleration profile of a minimum time retar- 
geting maneuver. The commanded angular acceleration profile was: 

. / v f 0.3 rad.s" 2 if 0 < t < .5(s); 

U ' \ -0.3 rad.s~ 2 if.5(s) < t < 1(5); 

Presumed angular accelerometer and gyro noise covariances were trans- 
formed into equivalent process noise for the feedback-linearized filter, with 
the additional simplification of neglecting a square noise term correspond- 
ing to the second entry U 2 = u> 2 of the equivalent input u. Presumed strain 
gauge sensor noises were taken from the literature. The sensor noise covari- 
ance matrix was modeled as a diagonal matrix with all diagonal elements 
equal to 0.00018 .Likewise the input covariance for U\ was 0.000005 and for 
U3 was 0.00001 . The input covariance related to the input U 2 was supposed 
to be negligible with respect to the other two. Two simulations were made 
, each one of them to illustrate the two methods of estimation described 
above.The first simulation was made without active bias suppression, i.e 
f = 0 in equation (10), the results being shown in Figures 3, 4. The sec- 
ond simulation was made with the use of bias suppression making b = 0 in 
eqaution (7) by letting f = -Kp,the results being shown in Figures 5, 6. 


4 Off Line Modeling 


In this section the method of Optimal Bilinear System Interpolation is used. 
In this technique the dynamical system is represented in in bilinear form (by 
active suppression ,if needed, of the bias term b), 


X =AX+B(X)u 
y = c 7 X 


where 2?(X) = [BiX | B2X | . . .]. This also means: 



• The (1/ 0) behavior is highly nonlinear. 


• The model is high dimensional if arising from Carleman linearization. 

• There is no clean ARMA model for system identification. 

Then, by using optimal interpolation one finds: 


• A closed form ,ci rcui t-implementable , reduced order, decoupled model 
which is also bilinear:cf. Figure 7. 


• (V 0)-based system identification can be used to tune the model to 
known responses to designer-selected typical excitations. 


• The dimension of the new system model is equal to the number of test 
signals. 


In the present application, the model is ”a priori” bilinear by the choice made 
for the inputs, so that the dimension is that of the structural model, here 
given by the number of mass points. For more realistic structural models, 
the filter dimension would nevertheless be high. Rather than tolerate the 
time delay found in the previous techniques of estimation, the method of 
operator spline interpolation can be used, to find the deflection amount 
between observations. The input-output (I/0)operator V, 

V : u -+ y (25) 

from the excitation vector u to an output vector y (such as v given by 
equation (2)) is imbedded in a Hilbert space of (I/O)operators of candidate 
bilinear systems, equipped with a reproducing kernel, see equation (31) 

t 

K t (u,v ) = exp J u T (s)R~ 1 v(s)ds (26) 
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where the weight matrix R is determined by eigenvalues of A in equation 
(24), which in the context of the present application corresponds to bounds 
on the structural frequencies. An interpolator of the form 

V t (u) = (27) 

i 

is constructed, tuned so that the structural responses to preselected test 
inputs u x are recorded, and optimally interpolating at system level the re- 
sponses to other excitations in the signal class. 

If the recorded system responses y x to the test input u x are reliably known, the 
’tuned’ coefficients C{ are obtained by solving the matrix equation 

G(t) c(t) = y(t) (28) 

where c(t) = co/{c t },y(£) = coI{y x (t)} and G(t) = {K t (u x ,ui)} l j. A more 
complex matrix equation yields the c, for uncertain y x :c f.[7], [8], [9]. The 
optimization leading to the functional interpolator Vj is formulated as a 
minimization of the maximum distance between the interpolating operator 
and any candidate operator that matches the experimental input-output 
signals. If the system data are not accurate, a weighted minimization that 
does not require exact matching of system responses can also be used. This 
minimization is carried out in a Hilbert space of input-output operators 
equipped with a weighted ”Fock space” scalar product which is the Hilbert 
sum of the causal (i.e., with lower triangular domains of integration) L 2 
scalar products of the kernels of the Volterra series of the operators in ques- 
tion, for which K% is the reproducing kernel. The general method is discussed 
in [9], although causality was differently implemented there, since symmetric 
kernels were used. 

The Hilbert space structure for m inputs(here m = 3) is defined as followsdet 

V.j •••,<«) = c T exp{(t- t 1 )A}B il ... 

. ..exp{(t n -i - t n )A}B in exp(t n A)X(0) (29) 

where B(X) = [B^X) | B 2 {X) \ ...} and 

These are the Volterra kernels for Vjj(fi), . . . ,Ui n (/ n ) so long as triangu- 
lar integration is used as in equation (30). Then the inner product is given 
as shown , 

< Vt, Vf >= Ejj . . . £ tn r tl . . . r, n f f ... 

Jo Jo 
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fin 

J Q t*'n(^ ^1 ) • • • ) tn)hn,ii,...,i n (t,t\ , . . . , t n ) 

dt n ...dti (30) 

with designer-selected weights r, > 0 corresponding to R = diag{r ,}, which 
yields the reproducing property: 


< V t ,K t (u,.) >= Vj(u) ( 31 ) 

The Volterra series for a bilinear system will yield a bounded norm < Vi, V t > 
provided the weights rj are chosen so that 

Si-'illB.f < j (32) 

where N = dimX and n a H is a bound on the real parts of the eigenvalues 
of A y so that ||e*p(,4*)l| 2 < N.exp(-2at) . The bound (32) is obtained af- 
ter cancellation of intermediate exponents in the factors exp{—2a(t{-\ 

which can be interchanged when computing L 2 bounds of |h n|ll ln (Mi, • • . ,t n )|, 

to guarantee negativity of the remaining exponential coefficients. The same 
bound is sufficient for having \y(t)\ = |V«(u)| -► 0 as t -> oo when u is 
L 2 (Square integrable) , as is found by use of the reproducing property (31) 
(Lower bounds are also needed when the inputs are not L ^ ,but are bounded 
almost everywhere:Dwyer,[7]). The advantages of such modeling are: 


• The model dimension is equal to the number of test inputs. 


• The modeling error is distributed throughout the chosen input sig- 
nal class (i.e by frequency or amplitude), rather than depending on 
nearness to a single reference input. 

• The interpolated signal(response) can be proven to converge asymp- 
totically to the true system response for any (unknown) excitation in 
the chosen signal class. 

In this technique of modeling, the real data f]^ can be recorded by ex- 
citing the real system with (constant or nonconstant )test inputs to construct 
the interpolator. The test inputs can be chosen to approximate the expected 
excitations of the system. Thus, the real system time responses are used for 
model matching, rather than responses synthesized from the mathematical 
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model. 

The problem with this technique, however, lies in the fact that storage of 
curves is required in order to compute the c,’s . The number of stored curves 
is equal to m x k x I N y where m is the number of test inputs, k is the di- 
mension of the output vector y, N the dimension of the state to be modeled 
and / is the number of possible initial values of each component. This diffi- 
culty does not allow the system to run in real time: e.g.,for the case of n 
point masses linked by massless but elastic connections one has N = 2n and 
l = n when measuring deflections or k = N — 2n if full state information is 
required in the planar motion case. 

5 Interpolator-Based Estimation: 

In this section, the two last techniques are combined to create a more effec- 
tive one by making use of the transition matrix spline of the bilinear system 
of the model: 

f)(t) = (33) 

In fact, the matrix- valued operator spline $ interpolates the transition ma- 
trices corresponding to the bilinear system model excited by constant 
or piecewise constant test inputs u^.This permits the construction of the 
response of the real time system in piecewise closed form, thereby replacing 
response curve storage by an analytic transition matrix generator, rather 
than the construction of the coefficients interpolator c t from the output test 
signals y T = (y 1 ,y 2 , . . . ,y m ). 

One gets each entry cf 9 of the matrix valued spline coefficients c,- by let- 
ting y' = in equation (28) where S^is the (p,qr)-entry in the transition 
matrix with constant u M : 

= exp{(A + J2 B 3 u< j ) ) t } ( 34 ) 

j 

The interpolated transition matrix is then used to update between obser- 
vations the structural state estimates obtained from an adaptive covariance 
filter based on a globally feedback-linearized transformation (seen in Section 
3.2) of the bilinear structural model. 

This last technique has the following features: 

• 4 is open loop, with fj made to match the real system at discrete 


intervals by re-initializing: 


($!>) =<t(M * ) (!(4)) 

(35) 

In contrast, the direct modeling of the I/O operator 


f) = ^ ~2,Ci(i)exp J u T R~ l u'dt 

(36) 


continuously tracks the true system time responses i^(t),but in this case 
c i(t) cannot be generated analytically and must be computed off-line. 

• The presence of an additive input does not give rise to a steady state 
tracking error observed in the earlier literature when additive as well as 
multiplicative inputs are present, as is the case for rapidly slewing struc- 
tures.Indeed, a convolution correction based on ♦ can be added, eliminating 
the need for active suppression of the bias term b in equation (7). 

• The number of curves to be generated is only m x N 2 instead of 
m x k x 7^, (where again N = 2n for the example of a structure composed 
of point masses connected by elastic appendages and in plane motion). 

• The possibly high dimensional recursive filter can run at a slower sam- 
pling rate chosen to be consistent with on board CPU capabilities. 

5.1 Interpolation Example 

An interpolator was designed for the same two bodies beam expander model 
previously described: The interpolator was optimized for input vectors u T = 
( u i) u 2i u 3) of the form (constant, 0,0), (0, constant, 0) and (0,0, constant) 

, chosen with a positive constant during the first half of a 1 second nominal 
minimum time rotation, and negative during the second half, for the first 
test input vector.The same positive constant was chosen throughout the 1 
second repointing for the second and third test input vectors, this qualita- 
tively correspond to the nominal angular motion where tl> i6 a square wave 
beginning at +0.3 and reversing, which yields positive (though not constant) 
values for uj 2 and 2 u :cf.Eq.(5). The constants were selected for boundedness 
of the interpolator according to [7], [8]. For the sake of convergence of the 
interpolator to the actual output of the system, and in the case of applying 
test inputs which are not square integrable ,the matrix R described above 
in the reproducing kernel analytical form should have diagonal elements 
satisfying the following condition: 


a/ m --/(a/m)2-||B,||2ATA,- 2 < 



r,-AT||J3;|| 2 < 

a/m + vW™) 2 - (37) 

li A is N X 7V,||i?i|| 2 = Y2pq i,pq wliere 6,- iM is the ( p,q ) element of B„ then 
A,, which is the upper limit of the i-th test input component, must be chosen 
to satisfy the bounds below, 

0 < A, < a/mVN\\Bi\\ (38) 

where a is an upper bound on the real parts of the eigenvalues of the matrix 
A ,and m is the number of test inputs applied to the interpolator. Indeed, 
the bound (41) allows selecting positive r[s in the inequality (37), which 
itself occurs in obtaining bounds on y(t) = V t (u) from equation (34), to 
guarantee y(t) — ► 0 :cf.[6] 

In case the applied test inputs are square integrable, as will be chosen on 
the example, one needs only a simpler form of bounds for r, described as 
follows, which implies the inequality (32) if all r, are the same, 

0 < riN\\B t \\ 2 < 2 a/m (39) 

The data used to drive this example are very near to those of a space 
based laser model: 

J 0 = 20,556Jfcs.m 2 ,m° = 10,720*5, m 1 = 152*5, p 1 = (0,14.421m) r 


( 1642 

121.6 \ 

c =( 162 

120 \ 

V 121.6 

1355.8 7 

’ V120 

1357 ) 


The nondiagonal form of the above matrices K and C is due to the fact that 
the two links that hold the small mirror in the top can be regarded as an 
equivalent link with equivalent stiffness and damping matrices ,as shown in 

[13]. 

Now by applying those data to the system and by the choice of constant = .3 
in the test inputs, one finds 0 < r, < 1/12 where r, = r for all t, giving 
R — diag(ri) = r[/] mxm .Two scries of simulations were made. In each 
one of them two alternatives , namely, bias elimination ,i.e., b = 0 or no 
external tip forces ,i.e, f = 0 are considered. In the first series of simulations 
the interpolator was used over the entire minimum time 1 second maneuver 
and thereafter. The results for the case of lateral deflection of the appendage 
are shown in Figures 8,9. It is important to notice here that in case of 
bias presence, a numerical convolution product was used with the original 
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transition matrix $ as well as with the interpolated transition matrix 
*(M*) for the sake of error comparisons. ’’Original” in the plots means the 
numerical solution of Eq.(l) , ’’discrete” means the numerical computation 
of 0)x(0) + $*(*, .) * b . One last numerical diffuculty was observed, in 
addition to care needed in generating the contribution from b by convolution 
(at least in this unusual beam expander example). That was the singular- 
ity of Eq.(28) at t=0 as already discussed in [9], which caused numerical 
unreliability during the first 0.25 seconds of this motion. Thus ’’transient 
error” is consistent with the ’’adaptive” nature of the interpolator which 
must learn from the system resonse to alleviate this numerical diffuculty, 
in the other series the state estimated in the example of section 3.3 was 
used to re-initialize the interpolator each .5 s The results are shown in Fig- 
ures 10,11. The re-initialization from the estimated state at given moments 
yielded better results than when using the same interpolator throughout. 

6 Concluding Remarks 

Filter alghorithms combine the propagation of measurements ’’between” ob- 
servations with updating of measurements ’’across” observations. Such up- 
dating to account for new observations has been shown here to be obtainable 
from an estimator based on a globally feedback-linearized model of a non- 
linear process. 

In case nonlinear transformation of the observed part of the process 
state is required , it was shown that the associated matrix Riccati differen- 
tial equation for the propagation of the estimation error covariance needs 
to be updated only in its ’’driving” term, given by the process noise co- 
variance . In contrast, all the coefficients of the Ricatti equation for the 
corresponding extended Kalman filter must be updated, so that consider- 
able CPU time is saved by pre— linearization, although the filter dimension 
is the same. Dimension reduction between measurements is therefore still 
desirable, motivating the next part of this work, which is reviewed below. 

If the process dynamics is ’’parametrically excited”, e.g. by gyroscopic cou- 
pling, it was then shown how the process state can be propagated between 
observations by interpolation of the input-output operator that maps the 
process excitation to the (time-varying) state transition matrix. This inter- 
polator was shown to be simultaneously optimized to match the measured 
system response to a set of pre-selected ’’test inputs”, which if chosen piece- 
wise constant can also be encoded analytically in closed form. 

Since the interpolator dimension is determined by the dimension of the pro- 



cess state ,it is therefore faster to periodically re-start the interpolator at 
the rate the feedback linearized filter( or any other) can process full state in- 
formation, rather than be tied down to the full order filter processing rate. 
New results were then given on the application of the method to the on- 
line estimation of transverse deflections of a rapidly slewing, gyroscopically 
slew-coupled beam expander, previously reported in [12]. Applications to 
line of sight disturbance error bound estimation in sliding control following 
[10] are also outlined in [11]. the interpolation technique may also be used 
to update observations for multi-axis motion of multibody systems , but 
this latter work is still undeveloped. 
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Figure 2. 
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